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Abstract 

We derive and evaluate one- loop functional flow equations for the effective interactions, self- 
energy and gap function in spin-singlet superfluids. The flow is generated by a fermionic frequency 
cutoff, which is supplemented by an external pairing field to treat divergencies associated with 
the Goldstone boson. To parametrize the singular momentum and frequency dependences of the 
effective interactions, the Nambu interaction vertex is decomposed in charge, magnetic, and nor- 
mal and anomalous pairing channels. The one-loop flow solves reduced (mean-field) models for 
superfluidity exactly, and captures also important fluctuation effects. The Ward identity from 
charge conservation is generally violated, but can be enforced by projecting the flow. Applying the 
general formalism to the two-dimensional attractive Hubbard model, we obtain detailed results on 
the momentum and frequency dependences of the effective interactions for weak and moderate bare 
interactions. The gap is reduced by fluctuations, with a stronger reduction at weaker interactions, 
as expected. 
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I. INTRODUCTION 



Numerous interacting Fermi systems undergo a phase transition associated with spon- 
taneous symmetry breaking at sufficiently low temperatures. Mean-field theory captures 
salient features of the symmetry-broken phase such as long-range order, quasi-particle exci- 
tations, and collective modes. For example, the BCS wave function provides a surprisingly 
faithful qualitative description of the superfluid ground state of an attractively interacting 
Fermi gas not only at weak, but also at strong coupling.^ Nevertheless, fluctuations often 
play an important role, both above and below the energy scale for symmetry breaking. At 
high energies, they renormalize the effective interactions generating an instability of the nor- 
mal (symmetric) state, which may enhance or reduce the scale for symmetry breaking. At 
low energies, order parameter fluctuations usually suppress the order at least partially. Trig- 
gered by the possibility of designing tunable attractively interacting Fermi systems in cold 
atom traps, the issue of fluctuation effects in fermionic superfluids has attracted renewed 
interest P 

A framework to deal with fluctuation effects on all energy scales is provided by the 
functional renormalization group (fRG). This method provides a flexible source of new ap- 
proximation schemes for interacting Fermi systems,^ which are obtained by truncating the 
exact functional flow for the effective action as a function of a decreasing infrared cutoff 
\\MI\ The common types of spontaneous symmetry breaking such as superconductivity or 
magnetic order are associated with a divergence of the effective two-particle interaction at a 
finite scale A c in a certain momentum channelP^ To continue the flow below the scale A c , 
an appropriate order parameter has to be introduced. 

A natural procedure is to decouple the interaction by a bosonic order parameter field, via 
a Hubbard-Stratonovich transformation, and to study the coupled flow of the fermionic and 
order parameter fields. Thereby order parameter fluctuations and also their interactions can 
be treated rather easily. This route to symmetry breaking in the fRG framework has been ex- 
plored already in several works on antiferromagnetismpEI' anc i superconductivity.^^ Typ- 
ically the bare microscopic interaction can be decoupled in various channels. An effective 
interaction with only one bosonic field corresponds to a strongly simplified representation 
of the effective two-fermion interaction. Systems with competing instabilities corresponding 
to distinct order parameters require the introduction of several bosonic fields.^ESl 
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Alternatively one may explore purely fermionic flows in the symmetry-broken phase. 
This can be done by adding an infinitesimal symmetry breaking term to the bare action, 
which is promoted to a finite order parameter below the scale A c . 2 * A simple one-loop 
truncation of the exact fRG flow equation with self-energy feedback was shown to yield 
an exact description of symmetry breaking for mean-field models such as the reduced BCS 
model, although the effective two-particle interactions diverge at the critical scale A C J 21 * 22 I A 
subsequent application to the two-dimensional attractive Hubbard model showed that the 
same truncation, with a rather naive parametrization of the effective two-particle vertex, 
yields surprisingly accurate results for the superconducting gap at weak coupling. 23 However, 
the flow could be carried out down to A = only for a symmetry-breaking pairing field Ao 
above a certain minimal value. At that value a spurious divergence of the two-particle 
vertex was found. Fortunately, the minimal Ao was rather small, more than two orders of 
magnitude below the size of the gap at the end of the flow, and in this sense close to the 
ideal case of an infinitesimal A . 

In this paper we further develop the fermionic fRG for spin-singlet superfluids as a pro- 
totype for a broken continuous symmetry. We stay with the one-loop truncation used pre- 
viously, but we derive and apply a much more accurate parametrization of the momentum 
and freqency dependence of the flowing two-particle vertex, taking all singularities in the 
particle-particle and particle-hole channel into account. We build on recent work on the 
structure of the Nambu two-particle vertex in a singlet superfluid,^ where constraints from 
symmetries (especially spin-rotation invariance) were derived, and insight into the singulari- 
ties associated with superfluidity was gained by analyzing the exact fRG flow of a mean-field 
model with charge and spin forward scattering in addition to the reduced BCS interaction. 
Furthermore, a decomposition of the Nambu vertex in distinct interaction channels was de- 
rived, extending the decomposition formulated by Husemann and Salmhofer^for the normal 
stated which will now be used to separate regular from singular momentum and frequency 
dependences. With an adequate parametrization of the vertex at hand we can fully explore 
the performance of the one-loop truncated fermionic RG for symmetry-breaking beyond 
mean-field models. Explicit results for the effective interactions, the self-energy, and the 
gap function will be presented for the two-dimensional Hubbard model with an attractive 
interaction as a prototypical case. The Ward identity relating the gap to the vertex in the 
phase fluctuation channel (Goldstone mode) is not consistent with the truncated flow. The 
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deviations are small at weak coupling, but they increase with the interaction strength. This 
problem can be treated by projecting the flow on the manifold of effective actions which 
respect the constraint imposed by the Ward identity We also analyze to what extent effects 
of the Goldstone mode on other channels are captured by the one-loop truncation. 

The article is structured as follows. In Sec. II the basic one-loop flow equations for the 
self-energy and the Nambu two-particle vertex are written down. Symmetry properties of the 
Nambu vertex following from spin rotation invariance and discrete symmetries are reviewed 
in Sec. III. The channel decomposition for spin-singlet superfluids is derived in Sec. IV, 
and the general structure of the flow equations is discussed. In Sec. V the random phase 
approximation is revisited in the framework of the channel decomposed flow equations. The 
general formalism is applied to the attractive Hubbard model in Sec. VI, with results for 
the self-energy, the gap function and the effective interactions in all channels. Merits and 
shortcomings of the channel decomposed one-loop flow equations are summarized in the 
conclusions, Sec. VII. 

II. TRUNCATED FLOW EQUATIONS 

We analyze the superfluid ground state of attractively interacting spin-| fermions. The 
systems are specified by a bare action of the form 

Sty,$\ = - Yl[ih>-m\Mk« + Uty,$\ , (1) 

k,cr 

where ipka and ipka are Grassmann variables associated with creation and annihilation op- 
erators, respectively. The variable k = (k , k) contains the Matsubara frequency k in 
addition to the momentum k, and a denotes the spin orientation. £(k) = e(k) — ji is the 
single-particle energy relative to the chemical potential, and U[ip, i/>] describes a spin-rotation 
invariant two-particle interaction 

U[rj), V>] = \^[U (ki, k 2 , k 3 , h)5 aia4 5 a2U3 -U(k u k 2 , h, k 3 )5 aia3 5 a2<74 ] 

X ^fc 1 <r 1 ^fc a o a V'*3«sV'fc4<r4 ■ ( 2 ) 

Here and in the following the usual temperature and volume factors are incorporated in the 
summation symbols. 
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Our analysis is based on a truncation of the exact flow equation^ for the effective action 
r A [-0, -0], that is, the generating functional for one-particle irreducible vertex functions in the 
presence of an infrared cutoff A. The cutoff is implemented by endowing the bare propagator 
with a regulator function. r A [-0, -0] interpolates smoothly between the regularized bare action 
at the initial scale Ao and the final effective action V[ip, -0] in the limit A — > 0. Spontaneous 
breaking of the U(l) charge symmetry in the superfluid state can be treated by adding a 
small (ultimately infinitesimal) symmetry breaking field 

5<S[0, 0] = ^ [AoW^-k^frt + A*(&)0 fet 0_ fe J (3) 

k 

to the bare action, which is then promoted to a finite order parameter in the course of the 
flowP 

Expanding the exact functional flow equation for r A in powers of ip and 0, one 
obtains a hierarchy of flow equations for the n-particle vertex functions.^ We truncate the 
hierarchy at the two-particle level, including however self-energy corrections due to con- 
tractions of three-particle terms.^ This truncation was used in all previous fermionic fRG 
studies of symmetry breakingP^^l j s eX act for mean-field models. 

In a superfluid state it is useful to use Nambu spinors (f>ks and (fiks defined as 

(f) k+ = -0fet, 4>k+ = fpkt, 4>k- = tp-hi, 4>k- = ^-fei (4) 

instead of ipka and ipka as a basis. The effective action as a functional of the Nambu fields, 
truncated beyond two-particle terms, has the form 



k.S2 



k s 



1,S2 



S2 ( Pk3S3 ( Pk4S4 1 (5) 

fcl,...,fe4 S1,...,S4 

where r^°^ A yields the grand canonical potential.^ For spin-singlet pairing with unbroken 
spin-rotation invariance only terms with an equal number of and fields contribute. 
The Nambu vertex r^i A S3S4 (A;i, fc 2 , k 3 , k^) is non-zero only for k% + k 2 = k 3 + k A , due to 
translation invariance. The (scale-dependent) Nambu propagator G A is related to r^ 2 ^ A by 
(G A ) _1 = r^ 2 ^ A , and can be written as a 2 x 2 matrix 

G * (fc)= |GU*><J}-wW^w ^ (6) 

G\{k) G A _(k) / I F* A (k) —G A (—k) ' 
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FIG. 1: Diagrammatic representation of the flow equation for the self-energy. The slashed line 
represents the single-scale propagator S A . 

The anomalous propag ator F A (k) is a symmetric function of ko and k. The Nambu self- 
energy E A is defined by the Dyson equation (G A ) _1 = (G A ) _1 — E A , where G A is the bare 
regularized propagator (in presence of A ). In the superfluid state it has the form 

I ^ A (k) A (&) — A A (k) \ 
S A (fc) = K } ° { J \ , (7) 

I Aq(/c) — A* A (&;) -£ A (-£;) I 

where S A (fc) is the normal component of the self-energy and A A (/c) is the (flowing) gap 
function. 

The gap function and the Nambu vertex are related by a Ward identity following from 
global charge conservation (see, for example, Ref. [21]) 



A\k) - Ao(fc) = [Ao(k')G A + (k')G A s ,(k') — Al(k')G^_(k')G J l a ,(k') 

k' s.s' 

xrf s t(k,k>,k',k). (8) 



The Ward identity implies that some components of the Nambu vertex diverge in case of 
spontaneous symmetry breaking (A A finite for A —> 0), which is a manifestation of the 
massless Goldstone boson. 

The flow equation for the Nambu self-energy, shown graphically in Fig. 1, is given by 

^EE s «w r i(^'^) > ( 9 ) 



where 



S A fixed 



[1 - G A (A;)£ A (A;)] 1 ^® [l - S A (A;)G A (A;)] 1 (10) 
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FIG. 2: Diagrammatic representation of the flow equation for the two-particle vertex. The dots 
denote differentiation of the propagator products with respect to the scale A. 

is the so-called single-scale propagator. The truncated flow equation for the Nambu vertex 
(see Fig. 2) reads 

Arg^fa, k 2 , k 3 , h) = nJSUCh, k 2 , k 3 , h) - nJS&jk, k 2 , k 3 , h) 

2 ' "SI S :;-M ' 



-TF s *{kiMMM) , (ii) 



where 



dA 1 

x r?$^fap,QMr$£^{qMM , (12) 



dA 1 

x r^ S4 (^,p,?^4)rg ) s f S34 ( g ,A ;i ,A;3,p) , (13) 
n^Jfci.k,**,**) = EE ^[^(p)^(9)] 

p,q si,...,a^ 

x fc, g,p)r^ S3S4 (p, g, fcs, fc 4 ) . (14) 

The flow equation for the self-energy is exact (for an exact r^ 4 ^ A ), while in the flow of 
p( 4 ) A contributions from r^ 6 ^ A beyond self-energy feedback have been discardedP^Zl These 
discarded contributions are at least of order (r^ A ) 3 , and they involve overlapping loops 
leading to a reduced momentum integration volume. The truncation is exact for mean-field 
models with a reduced BCS and/or forward scattering interaction, although r^ A becomes 
large at the critical scaleP^^ The particle-particle terms in the Nambu representation 
contain particle-hole contributions in the original fermion basis and vice versa. In particular, 
the usual particle-particle contribution driving the Cooper instability is contained in the 
Nambu particle-hole diagrams. 



III. SYMMETRIES OF NAMBU VERTEX 



The Nambu vertex T^J^ S 3s 4 (ki, k 2 , k 3 , k 4 ) has 16 components corresponding to the choices 
Si — ± for i = 1, ... ,4. Spin rotation invariance reduces the number of independent com- 
ponents of the Nambu vertex substantially. In Ref. G31 it was shown that the vertex can 
be parametrized by three functions of k\, k 2 , k 3 , k 4 , where k\ + k 2 = k 3 + k 4 for transla- 
tion invariant systems. These functions are further constrained by discrete symmetries. In 
this section we describe the spin-rotation invariant form of the Nambu vertex as derived in 
Ref. EH 

In addition to the normal interaction involving two creation and two annihilation op- 
erators (ipipi/jip) there are also anomalous interactions corresponding to operator products 
+ conjugate and ^ipipip + conjugate. 1 21123 1 We first write down manifestly spin-rotation 
invariant forms for each of these terms in the 7/>-basis, and then translate to the Nambu rep- 
resentation. 

A spin-rotation invariant normal interaction can be written as^ 



r (2+2) [?/>, ip] = ~y~] [V(h, k 2 , k 3 , k i )5 aia4 5 a2(T3 - V(kt, k 2 , k 4 , k 3 )5 aiCT3 S a2(T4 ] 

ki ,(T^ 

X i'k 1 a 1 4'k 2 a 2 ^k 3 a 3 ^k i a 4 ■ (15) 



In this section we suppress the superscript A denoting the scale dependence. One may also 
write r( 2+2 )[^,^] 

as a sum of a spin singlet and a spin triplet components 



T {2+2) [ip, ip] = - ^2 [ vS (ki, h, k 3 , h) S ai(T2a3<T4 + V T (k 1 , k 2 , k 3 , fc 4 ) T aicr2<T3a4 \ 

ki ,(T i 

x i)k^3k^h z ^k^ A , (16) 

V s (k h k 2 , k 3 , fc 4 ) = V(k 1: k 2 , k 3 , k 4 ) + V(h, k 2 , h, k 3 ) , (17) 
V T (k u k 2 , h, k 4 ) = V(k u k 2 , k 3 , h) - V(k u k 2 , fc 4 , k 3 ) . (18) 



Time reversal invariance^ and Osterwalder-Schrader positivity (corresponding to hermiticity 
of the underlying Hamiltionian) yield the following relations: 

V{k u k 2 ,k 3 ,k A ) = V(Tlh,'Jlk 3 ,'Jlk 2 ,'Jlk 1 ) , (19) 
V(ki, k 2 , k 3 , k 4 ) = V*(kl, k* 3 , k* 2 , kl) , (20) 
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where TZk = (k , — k) and k* = (— k ,k). Spatial inversions also transform k to TZk, but 
without permuting the momenta in V(k±, k 2 , k 3 , k 4 ). The invariance of the action under 
the exchange of identical particles yields separate symmetry relations for V s and V T under 
k\ -H- k 2 and k 3 -H- fc 4 : 

^(^,^,^3^4)= V s (k 2 ,h,k 3 ,k 4 ) = V S (h, k 2 , k 4 , k 3 ) , (21) 
^(fci, fc 2 , A; 3 , fc 4 ) = -^ T (A; 2 , k u k 3 , k 4 ) = -^(fci, fc 2 , *&) , (22) 

while V obeys only 

^(fci, k 2 , h, k 4 ) = V(k 2 , k u k 4 , k 3 ) . (23) 

A spin-rotation invariant anomalous interaction involving four creation or four annihila- 
tion operators can be written in the form 

r( 4 +°)[^,^] = 

Id 

- w T (h, fc 2 , k 3 , h) [(4>krf4>k 2 i + i>k 1 ii>k 2 ^)(i>k^i>k 4 i + 4>k 3 i4>k^) 

-2(^fc 1 t^fc 2 t^*34.^fc44- + ^k^k^k^k^)] + conj.j , (24) 

consisting of a singlet (S) and a triplet (T) term, which are separately invariant under spin 
rotations. Conjugated terms denoted by "conj." are obtained by reversing the order of 
fields, replacing ip^ by ipk*a, and complex conjugation of the functions W ,T . Time reversal 
invariance yields a relation between W S,T and W S ' T *, which assumes the simple form 

W s,T *(k 1 , k 2 , k 3 , k 4 ) = W s ' T (-k u -k 2 , -k 3 , -k 4 ) , (25) 

if the order parameter (gap function) is chosen real. Furthermore, the functions 
W s (ki, k 2 , k 3 , k 4 ) and W T (ki, k 2 , k 3 , k 4 ) are symmetric and antisymmetric under the ex- 
change k\ -H- k 2 or k 3 -B- k 4 , respectively, and symmetric under the simultaneous exchange 
ki -H- k 3 and k 2 -H- k 4 : 

W s (k u k 2 ,k 3 ,k A ) = W s (k 2 ,k u k 3 ,k 4 ) = W s (k u k 2 ,k 4 ,k 3 ) 

= W s (k 3 ,k 4 ,h,k 2 ) , (26) 

W T (k 1 , k 2 , k 3 , k 4 ) = -W T (k 2 , h, k 3 , k 4 ) = -W^, fc 2 , fc 4 , fc 3 ) 

= W T (k 3 ,k 4 ,h,k 2 ) . (27) 
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Finally, a spin-rotation invariant anomalous interaction with three creation and one an- 
nihilation operators, or vice versa, can be written as 

hi cr 

X T (k 1 , k 2 , k 3 , fc 4 ) J ^ ^k^i^^ksi + ^k 2 l^k 3 t)^k 4 a 



+ 



+ conj 



}■ 



(28) 



where e-f = 1 and e^. = —1. The terms with the coefficients X 5 and X T are separately 
spin-rotation invariant. Time reversal invariance yields the relation 

X s ' T *(k 1 , k 2 , k 3 , k 4 ) — X s ' T (—k 1 , —k 2 , -k 3 , -k 4 ) (29) 

for a real choice of the order parameter. Invariance under particle exchange yields 



X s (k u k 2) k 3 , k 4 ) = X s \k u k 3 , k 2 , k 4 ) , 
X T (ki,k 2 ,k 3 ,k 4 ) = -X T (ki, k 3 , k 2 , k 4 ) . 



(30) 



To express the two-particle interaction in terms of Nambu fields, it is convenient to collect 
the 16 components of the Nambu vertex rii S2S3S4 in a 4 x 4 matrix 



r (4) 



r +++ 


F (4) 
1 ++-4 


F (4) 
- 1 +-++ 


F (4) 

1 + + 


r (4) 

1 +++- 


F (4) 
1 ++ — 


F (4) 


F (4) 
1 + 


r (4) 

1 -+++ 


F (4) 
1 -+-4 


F (4) 

- 1 ++ 


F (4) 

1 + 


r (4) 

1 -++- 


F (4) 
1 -+ — 


F (4) 

1 +- 


p(4) 



(31) 



/ 



Note that the rows in this matrix are labeled by s\ and s 4 , while columns are labeled by s 2 
and S3. With this assignment the Bethe-Salpeter equation yielding the exact Nambu vertex 
in reduced (mean-field) models can be written as a matrix equation. Translating the spin- 
rotation invariant structure of the various interaction terms described above to the Nambu 
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representation, one obtains a Nambu vertex of the following formpD 
T^(k 1 ,k 2 ,k 3 ,k 4 ) = 

' V T (k u k 2 ,k 3l k 4 ) X(k u k 2 ,k 3 ,k 4 ) X*(kt,ktkt,kt) -V(k u -k 3 ,-k 2 ,k 4 )^ 
-X(ki,k 2 ,h,k 3 ) W(k!,k 2 ,k 3 ,k4) V(fa, -k 4 , -k 2 , k 3 ) X*(k 4 ,k 3 ,k 1 ,k 2 ) 
-X*(kt, k* 3 , k 2 , k\) V*{k x , -h, -k 2 , k 3 ) W*(kt k* 3 , k 2 , fcf) X(kl, k*, kt k* 4 ) 

y—V*(ki, —k 3 , —k 2 , k 4 ) -X*(k 4 , k 3 , k 2 , h) -X(kl,k*,k*,k$) V T * (k u k 2 , k 3 , k 4 ) j 

(32) 

The matrix elements W and X are related to the anomalous (4+0) and (3+1) interactions, 
respectively: 

W(h,k 2 ,k 3 ,h) = W^k^-h.-k^k^-W^ku-k^-h,^) 
+ W T (k 1 , -k 4 , -k 3 , k 2 ) - W T (k 1 , -k 3 , -k 4 , k 2 ) 
+ 2W T (k u k 2 ,-k 3 ,-k 4 ) , (33) 

X(k!,k 2 , k 3 , k A ) = X s (k l ,k 2 ,-k 3 ,k 4 ) - X s (k 2 ,k 1 ,-k 3 ,k 4 ) 
+ X T (k 1 , k 2 , -k 3 , k 4 ) - X T (k 2 , fci, -k 3 , k 4 ) 
+ 2X T (-k 3 ,k 2 ,k u k 4 ) . (34) 
The functions W and X obey the following relations under exchange of variables: 
W(k u k 2 , k 3 , k 4 ) = -W(k 2 , k u k 3 , k 4 ) = -W(k u k 2 , k 4 , k 3 ) 

= W(-k 4 ,-k 3 ,-k 2 ,-h) , (35) 

X(k u k 2 , k 3 , k 4 ) = -X(k 2 ,ki, k 3 ,k 4 ) . (36) 

For translation invariant systems the functions V(ki, k 2 , k 3 , k 4 ), W(ki,k 2 ,k 3 ,k 4 ) and 
X(k 4 , k 2 , k 3 , k 4 ) are non-zero only if k\ + k 2 = k 3 + k 4 , and can therefore be parametrized by 
three energy and momentum variables. Reflection invariance implies that these functions 
are invariant under kj y — kj. 

IV. CHANNEL DECOMPOSITION 

The two-particle vertex acquires a pronounced momentum and frequency dependence in 
the course of the flow, which becomes even singular at the critical scale for spontaneous 
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symmetry breaking. A parametrization based on weak coupling power counting is not ad- 
equate in this situation. Keeping the full dependence on the three independent momenta 
and frequencies is technically not feasible. The particle-particle and particle-hole contribu- 



tions to the flow, Eq. (11), depend strongly on certain linear combinations of momenta and 
frequencies, namely 



n PH,d 



n 



PH,cr 



n 



pp 

S1S2S3S4 



(k 1 ,k 2 ,k 3 ,k 4 ) 
{kx,k2,k 3 ,k4) 
(^1,^2,^3,^4) 



k 3 -k 2 , 
h - h , 
h + k 2 ■ 



(37) 



This is because the poles of the contributing propagators coalesce when the above combi- 
nations of momenta and frequencies vanish or are situated at special nesting points (in case 
of nested Fermi surfaces). We therefore write the vertex as a sum of interaction channels, 
where each channel carries one potentially singular momentum dependence, which can be 
parametrized accurately, while the dependence on the remaining two momentum variables 
is treated more crudely. This channel decomposition was introduced by Husemann and 
Salmhofer^ for the two-particle vertex in a normal metallic stated and extended by us for 
a superfluid stated Most recently it was also formulated for an antiferromagnetic stated 



A. Interaction channels 

Following Husemann and Salmhofer,^ we write the normal vertex in the form 

r( 2+2 )>,v>] = 

1 

2 



+ 7T C i t 1 +k A k 2 +k 3 {k 3 — k 2 ) ^kxa^k^a'^k-iu'^k^a 

t 2 ' 2 / 

Vj M ^ 1+fc4 k 2 +k. i (k 3 — k 2 ) / JT(j 1(J4 • To- 2(T3 V'fci 0-1^20-2 ^30-3 ^4 c 



2 

ki CTi 

+ „ P i\-k2 k 4 -k :i {h + k 2 ) V" 4> kl a4 ! k 2 a>4 ! k 3 a>4 ! k i a , (38) 

1 u 2 ' 2 / 

ki <T,cr 

where U[ip, ip] is the bare interaction, and the coupling functions C A , M A , and P A capture the 
"charge", "magnetic" (spin), and "pairing" channels, respectively. The matrices collected 
in r = (t x ,t v ,t z ) are the three Pauli matrices. The prime at the sums over ki indicates 
momentum (and frequency) conservation, ki + k 2 = k 3 + k 4 . The momentum argument 
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in brackets is the momentum transfer for the charge and magnetic channels, and the total 
momentum for the pairing channel. These are the variables for which a singular dependence 



is expected. Comparing the ansatz Eq. (38) to the general spin-rotation invariant form of 
the normal vertex Eq. (15), written in terms of V A (k\, k 2 , k 3 , k±), one obtains the relation 

V A (k 1 ,k 2 ,k 3 ,k 4 ) = U(k lt k 2 ,k 3 , k 4 ) 



+ 



C k 1 +k 4 k 2 +k 3 (k 3 — k 2 ) + P k 1 -k 2 k 4 -k 3 {k\ + k 2 ) 



M 



k ±+k 4 k 2 +k 3 
2 ' 2 

/1/fA n^rl PA 



2M kl +k 3 k 2 +k 4 (k\ — k 3 



ki+k 2 ,k 3 +k 4 



■ (39) 



The flow equations for C , M and P are obtained by choosing a Nambu component 
involving the normal interaction V A , such as F+_ + _(ki, k 2 , k 3 , /c 4 ) = V A (ki, — fc 4 , — k 2 , k 3 ), 
and linking the flow of the various components to the Nambu particle-particle and particle- 
hole terms such that momenta in brackets correspond to the strong momentum dependences 



as in Eq. (37). One thus obtains^ 

1 



d_ 
dk 

d 

d 



MtAq) 



4 



-n 



pp 



n 



PH,cr 

+-+- 



\k -\- 2 , 2 kj k -\- 2 , 2 ^ ) 

(h 4- i —l — P k — 2 2 

\^ ~ 2' 2 ' 2 ' 2 



k') 



n 



pp 



-|- o j n k^ k -\- 21 'o k 



Pkk'io) — n + _' + _(A: + § > & 2'^+2'^ I) ■ 



(40) 
(41) 

(42) 

p&{q) = Pffb) + P&b) , (43) 

where P fc ^, (g) is symmetric under sign changes of k and A;', while P kk \ (g) is antisymmetric. 
For the anomalous (4+0)-interactions the dependence on the (total) momentum of the 



^Y'/c/s-Vi/ "1 1 — V" 1 2' 2' 1 2' 2 

The pairing interaction can be split into a singlet and a triplet component as 

jA /„\ _ n>SAf~\ , oT,A 



Cooper pairs contained in T^ +0 ^ A [ip,ip], Eq. (24), is expected to become singular, which is 
taken into account by the ansatz 



W s > A (h,k 2 ,k 3 ,h) = Wt[ 



S,A 



-k 2 k 4 -k 3 



k3 (h + k 2 ) S kl 



+k 2 +k 3 +k 4 ,0 



W T ' A (ki, k 2 , k 3 , h) = W j;\ ki _ k3 (fci + k 2 ) 5 kl+k2+ka+k4fi . 



2 

T T,A 



(44) 
(45) 



Eq. (33) then yields 



W A (h,k 2 ,k 3 ,h) 



+ w 



r S,A 

+1 

2 

T,A 



W^l k4 k 2 +k 3 (k 3 — k 2 ) 



^ fci'+fc 3 fc 2 +- i - h 1 

2 > 2 



; i I i i ,v , . ;. . ( k 3 k 2 ) I 1 . ; / j | /. , 



rT,A 

2 



(h - k 3) 
{h - k 3 ) 



k 2 k 3 -k 4 
> 2 



(fci + &2 



ki+k 2 ,k :j +k 4 



(46) 



13 



A Nambu vertex component capturing this interaction is T+^__(fa, fa, fa, k±). Matching 
again the strong momentum dependences in brackets with those of the particle-particle and 
particle-hole terms, one gets 24 



_d_ 
dA 



r PH,d 



dA 



__n pp 



.(& + §>! 



k, I k j k ~\~ 



21 ' 



3 £ _ h & — h* h ! + 2" 

n,, o rv j r\j \ 2 J 



(47) 
(48) 



The functions X s ' A (ki,fa,fa,fa) and X T,A (fci, &2, &3, £4) parametrizing r^ 3+1 ^ A [^, ■?/>] in 



Eq. (28) are expected to depend singularly on &2 + &3, which is the total momentum of the 



Cooper pair contained in T^ 3+1 ^ A [ip, ip\. We therefore write 



X s > A (fa,fa,fa,fa) 
X T,A (fa, k 2 , fa, fa) 



X 
X 



(fa 



S,A 

fc l+ fc 4 k 2~ k 3 
2 ' 2 

T,A /, 

fcX+fc 4 k 2 -k 3 \™2 
2 ' 2 



^3) 3k 1 +k2+k 3 ,k4 j 
^3) ^fcl+fc 2 +fc3,fc 4 • 



(49) 
(50) 



Eq. (34) then yields 



X A (fa,fa,fa,fa) — X f[+ fc4 fc2+fc3 (A:2 — fa) — -Af f A fc4 ki+kg (fa — fa 

2 ' 2 2 ' 2 

+ X k[+k A k 2 + k 3 (fa ~ fa) ~ X k ' 2+kA fcl+fc3 (fa — fa) 

2 

5 



2 ' 2 



+ fc2 _ fc] (fc! + fc 2 ) 



fcl+fc2,fc3+fc4 



(51) 



Anomalous (3+l)-interactions are contained in the Nambu vertex component 
T+]^_ + (fa, fa, fa, fa). Matching singular momentum dependences between the vertex 
on the left hand side and the particle-particle and particle-hole terms on the right hand 
side of the flow equation yields^ 



d 

dA 



PH,d 



2' ^ ~ 2i 



n_i__i y(k ~ I - 2 ' 2 k > o k,k ~\~ ' 



■ 2 



2) ' 



d - r T,A 



k,k + 



^Y"fcfc' ^ 4"++— t-V" 1 2' 2 u '2 '">'" 1 2' 



(52) 
(53) 



Eqs. (40)- (42), (47), (48), (52), and (53) determine the flow of the complete set of coupling 

A A A S A. T j\ 

functions describing the Nambu vertex, that is, C kk ,(q), M kk ,(q), P kk ,(q), W kk , (q), W kk , (q), 



X kk , (q), and X kk f(q), respectively. Note that the above choice of Nambu components is not 
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unique. Any component containing V A , W A , and X A , respectively, could have been chosen. 
The resulting equations for the functions C kk ,(q) etc. are equivalent. 

Discrete symmetries and Osterwalder-Schrader positivity constrain the functions C kk , (q) 
etc. by relations analogous to those for the interaction functions presented in Sec. III. The 
normal interaction components obey 

cUq) = c A mk ,(Kq) = c A k ,{- q ) = c%_ k ,( q ) , (54) 

= KknAKq) = M A ,(-q) = M^-M , (55) 
pUq) = CM = P k>M = P-l',-k(-<l) , (56) 

where the first equation follows from inversion symmetry, the second from inversion and 
time reversal symmetry, and the third from inversion symmetry and positivity. For the 
anomalous (4+0)-interaction, invariance under spatial inversion and time reversal yield the 
relations 

Wtf(q) = W^iKq) = W£_ w {-q) (57) 
for v = S,T, and for the (3+l)-interactions 

K k %) = XZLsW) = Xt-A-Q) ■ (58) 
The complete Nambu vertex can be written in the form (see Fig. 3) 

(h, k 2 ,k 3 , k A ) 



+ 



K?£s< k ^h-h)- V^ A S4 (^, k 3 - ki) 

■\-^:t: h + h)} s kl+k2 , k3+k4 , (59) 



where the first term represents the Nambu components of the bare interaction, while the 
other terms are generated by the particle-hole and particle-particle contributions to the flow, 
that is 

^KZLA^, ^ k * - h) = n™> 2 % Si (h, h, h, k 4 ) , (60) 
^v^^^-M + h) = -\nZ 2S:iSi (kuk 2 ,k 3 ,h) . (6i) 

The crossed particle-hole contribution yields the flow of V FK,A with indices 1 and 2 exchanged 
and a minus sign compared to the direct contribution. Collecting terms with the variable 
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*4 \ 



m 



V 



X 




FIG. 3: Decomposition of the Nambu vertex in bare interaction, particle-hole channels, and 
particle-particle channel. 

q = &3 — ki in the channel decomposition, and writing the components in matrix form as in 

/ 



Eq. (31), one obtains 



V 



PH,A 



XLi-q) X£ k ,(q) 



-A 



W&(q) P£M 



(62) 



where 



XtAq) = X S J(q)+X^(q), 



and 



K^(Q)=Ch,(q)±M£ k ,(q). 
Collecting terms with the variable q = k\ + A; 2 , one finds 



(63) 
(64) 

(65) 



V 



PP,A 



(k,k';q) 



( p T J(q) 


~Xlt(q) 


~X^(q) 


KM 




-W^(q) 




~X T J k %(q) 


x T iM 




~W T J\q) 


-X T J k %{q) 




X^iq) 


xlt{ q ) 





\ 



(66) 



Note that V PH,A captures the full information on the coupling functions C kk ,(q) etc. char- 
acterizing the Nambu vertex. By contrast, V PP,A collects only magnetic and triplet pairing 
components. 

For q = the matrix V PH ' A has the same structure as the Nambu vertex for a mean-field 
model with reduced BCS and forward scattering interactions.^ Contributions with q ^ 
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correspond to fluctuations away from the zero momentum Cooper and forward scattering 
channels. 

It is convenient to use linear combinations of P A and W A corresponding to amplitude 
and phase variables. For a real gap function these combinations are 

a a M = Re[p A ,(q) + w A ,( q )} , (67) 

Qh(q) = Re[P A ,(q) - W A ,(q)} ■ (68) 

Amplitude and phase variables for singlet and triplet components can be defined by analo- 
gous linear combinations. Note that P kk , (q) and W kk , (q) are generally complex functions for 
q ^ 0, even for a real gap. For their real and imaginary parts we use the notation P kk /{q), 



W' kk ,(q) and P kk )(q), W kk )(q), respectively. Instead of the representation Eq. (31) it can be 
advantageous to use a Pauli matrix basis to represent the Nambu vertex, as described in 
Appendix A. 



B. Boson propagators and fermion-boson vertices 

To achieve an efficient parametrization of the momentum and frequency dependences, the 
coupling functions are written in the form of boson mediated interactions with bosonic prop- 
agators and fermion-boson vertices, as proposed by Husemann and SalmhoferP^ The bosonic 
propagators capture the (potentially) singular dependence on the transfer momentum and 
frequency while the fermion-boson vertices describe the more regular remaining momentum 
and frequency dependences. For example, the charge coupling function is decomposed as 

C£M = E C ^(9) 9L(k, q) 9 Uk', q) , (69) 

a,a' 

where the functions g A a (k,q) provide a real orthonormal basis set of fc-space functions, 
satisfying 

J dfx(k) g A a (k, q) g^,(k, q) = 5 aa , (70) 

with a suitable (not yet specified) fc-space measure dfi(k). Viewing C kk ,{q) as a boson 
mediated interaction, the functions C A a ,(q) can be interpreted as boson propagators and 
g A a (k,q) as fermion-boson vertices. Analogous decompositions are used for the magnetic 
and pairing coupling functions M kk ,(q) and P kk /(q), or the singlet /triplet components of the 
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latter. The anomalous (4+0) coupling function W kk , (q) can also be decomposed in the form 



Eq. (69). Alternatively one may decompose the amplitude and phase coupling functions. 
The anomalous (3+1) coupling functions X kk ,(q) require a more general decomposition 

X£M = J2 X -M9L(k,q)gL'(k',q) , (71) 

a,a' 

with two different sets of basis functions g£ a and g£ a , since the k and fc'-dependence of 
X kk ,(q) is generally different. 

Summing over a complete set of basis functions, the above decomposition is exact. In 
practice one has to approximate the infinite sum by a finite number of terms, with a suitable 
choice of boson-propagators and fermion-boson vertices. 



C. Classification of contributions to the flow 



Inserting the channel decomposed Nambu vertex on the right hand side of the flow 
equation yields several contributions which can be distinguished by their topology when 
representing the coupling functions by boson mediated interactions. For a graphical repre- 
sentation we use the symbolic notation from Fig. 3, where bosons mediating interactions in 
the (Nambu) particle-hole and particle-particle channels are represented by a wiggly and a 
double line, respectively. All contributions to the flow of the vertex are of second order in 
the interaction. We discuss the different topologies for diagrams with two wiggly lines as ex- 
amples. There are three distinct classes, which we refer to as "propagator renormalization" , 
"vertex correction", and "box contribution". For the propagator renormalization (Fig. 4, 
left) the momenta of both bosonic propagators coincide with the momentum transported 
through the fermionic bubble. Hence, a singularity in the bosonic propagators generated by 
the bubble is amplified by feedback of both propagators. For the vertex correction (Fig. 4, 
right) the momentum of one of the bosonic propagators coincides with the momentum of the 
fermionic (Nambu) particle-hole pair. Potential singularities of the other bosonic propagator 
are wiped out by the momentum integration. Note that at zero temperature all expected 
singularities of the vertex are integrable in two spatial dimensions, even the infrared singu- 
larity associated with the Goldstone mode. For the box contribution (Fig. 5) singularities 
of the fermionic pair are not amplified by singularities of the bosonic propagators which 
are both integrated. The contribution from the propagator renormalization diagram thus 
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FIG. 4: Examples for propagator renormalizaton (left) and vertex correction (right). The variable 
p is integrated. 




FIG. 5: Example for a box diagram with integration variable p. 

dominates in the formation of singularities at special wave vectors q = k 3 — k 2 . In mean-field 
models with reduced interactions it yields the complete flow, while vertex corrections and 
box contributions vanish. 

The assignment of momenta in the channel decomposition was designed to deal with 
singularities generated by the fermionic propagators. However, the box contribution exhibits 
another singularity generated by the singularity of the bosonic propagators at momentum 
zero. For k\ = k 3 (that is, k = k') the two bosonic propagators in Fig. 5 carry the same 
momentum variable. In the phase fluctuation channel (Goldstone mode) these propagators 
diverge quadratically at small momenta and frequencies (for A = 0). The product of 
two such singularities is not integrable in two dimensions. This problem can be treated by 
introducing a scale dependent pairing field Aq , which tends to zero continuously toward 
the end of the flow. A finite pairing field regularizes divergences in the Cooper channel 
(including the Goldstone mode), such that the right hand side of the flow equations remains 
finite at each finite scale, and the flow is integrable down to A — )■ 0, Aq — > 0, as discussed 
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A^^/l Ax^J'l 4^ » f 1 4^ » f 1 
2 **1> 3 3 <f^"V 3 >^^ii 2 2J ^ 13 

FIG. 6: Diagrammatic representation of contributions to the flow of y PH > A . The dot denotes a 
A-derivative acting on the product of the two fermionic propagators. 




1 4 

2 3. 




1 4 
2 3 





1 4 
2 3. 





FIG. 7: Diagrammatic representation of contributions to the flow of V 



PP,A 



in more detail in Sec. VI D 4 A scale dependent Ag does not modify the structure of the 



flow equations, it merely yields additional contributions to the scale derivative of the bare 
(Nambu) propagator G A . 

In addition to the contributions shown in Figs. 4 and 5 there are analogous contributions 
with wiggly lines replaced by double lines corresponding to the particle-particle channel 
and 4-point vertices representing the bare interaction (as in Fig. 3), including all possible 
mixtures of channels. The complete set of contributions to the flow of V PH,A and V FF,A is 
shown in Figs. 6 and 7, respectively. 
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V. RANDOM PHASE APPROXIMATION 



To gain insight into the singularity structure of the Nambu vertex it is instructive to 
consider the random phase approximation (RPA) before analyzing the full set of flow equa- 
tions. In the conventional formulation the RPA corresponds to a summation of all (direct) 
particle-hole ladder contributions to the Nambu vertex with bare interactions and mean- 
field propagators. The self-energy is obtained from the usual mean-field equation, that is, 
self-consistent first order perturbation theory. In the channel decomposed functional RG 
framework derived in Sec. IV, the RPA is equivalent to the approximation 

rg£ w (*, k'; q) = U slS2S3Si (k, k'; q) + V^Jk, k'; q) , (72) 

that is, the crossed particle-hole and particle-particle channels are discarded. Throughout 
this section we parametrize the momentum variables k\, &2, ks, k^ as ki — k+ |, k% — k' — §, 



&3 = k' + i, and ki = k—\. The flow equation (60 ) for V ' can then be formally integrated 



to obtain an integral equation which, expressed in term of r( 4 ) A , reads 
^s^s 3 s4,(k,k';q) = U slS2SaS4 (k,k';q) 



E E tW.«(*>J>; ?) G \4P ~ + I) rg^(p, k'; q) . (73) 



,(4)A 
P s't 

This is the familiar Bethe-Salpeter-type equation corresponding to a summation of (Nambu) 
particle-hole ladders. Using this equation, the flow equation for the self-energy Eq. ^ can 
also be integrated, yielding the usual mean-field equation 

Ks 2 = E E u & k '-> °) • ( 74 ) 

k' s' 1 ,s' 2 



The integral equation (73) can be written in matrix form such that Nambu index sums 
correspond to matrix products. In particular, choosing the Pauli matrix basis described in 
Appendix A, one obtains 

f( 4 ) A (A;, k'; q) = fj(k, k'; q) + E U(*,P5 ?) L A (p; q) f (4)A (p, k'] q) , (75) 

p 

where the components of L A (p; q) are given by 

l ii>te = \ E t&T& G t,(P - D G t 3 (P + I ) • (76) 
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For a spin-rotation invariant system, the bare interaction can be written in the form 

= ]- ^C^,{q)^ij k+q /2^k'- q /2,a4k'+ q /2,a4k- q /2,a 



2 

k,k',q 



2 

k,k',q Oi 

+ 2 P i°'^) 5Z^/ 2 + fe - CT ^/ 2 - fc ' CT '^/ 2 - fc, - CT '^/2+fc', CT > (77) 

k,k',q a,cr' 

which is analogous to the decomposition of the fluctuation contributions in Eq. ( |38| ). In the 
special case where the bare coupling functions C^,(q), M^,(q), and Pjfjiq) are non-zero 
only for g = 0, this becomes the reduced BCS and forward scattering interaction of the 
model discussed in detail in Ref. EH In that case the mean-field equation for the self-energy 
is exact, and the Bethe-Salpeter equation yields the exact vertex r^J^ 83S4 (k, k'; q = 0). 
For an explicit evaluation of the RPA vertex we assume separable interactions 

Cl%) = C^(q)f c (k)f c (k') , 
Ml%) = M^(q)Uk)Uk') , 

P { k%) = P { °\q)fp(k)f P (k') , (78) 

with symmetric (under k i— > —k) form factors, and a bare gap function A Q (k) = A f p (k) 
with the same form factor as the pairing interaction. The coupling functions contributing 



to V PH ' A , see Eq. (62), then also factorize: 

CUt) = C A (q)f c (k)fc(k') , 
M A k ,(q) = M k (q)Uk)f m (k') , 

pLM) = P A (q)f P (k)f P (k') , 

Wt«{q) = W A (q)f p (k)f p (k') , 

X A k ,(q) = X A (q)f c (k)f p (k>) , (79) 
and the gap function has the form 

A A (A;) = A A f p (k) . (80) 
The vertex assumes a particularly simple form in the Pauli matrix basis, namely 

t^ A (k,k>;q)={(k)t^ A (q){(k>), (81) 
22 



where i(k) is the diagonal matrix 



and 



with 



f (A:) = diag[/ m (A:), /„(*:), / P (A;), / C (A:)] 

f( 4 ) A (g)=U(g)+V PH ' A (g), 
/ 2M(°)(g) 



(82) 



(83) 



U(g) 



V 









and 



\™> A (q) 



( 2M A (q) 



V o 





p(0)( g ) 








\ 




2C^(q) J 



(84) 



A A {q) 



\ 

P" A (q) 2X' A (q) 

-P" A {q) $ A (g) -2X" A (q) 

2X lA {q) 2X" A (q) 2C A (q) J 



(85) 



Here A A (q) = P' A {q) + W' A {q) and $ A (g) = P' A {q) - W' A {q). Primes denote real parts 
and double primes imaginary parts. At go — all imaginary parts vanish. For q = the 
above matrix simplifies to the vertex previously obtained for the reduced BCS and forward 
scattering model, 24 in a slightly different basis yielding some sign changes. 



Inserting the factorized form of the vertex into the Bethe-Salpeter equation Eq. (75), one 
obtains a linear algebraic equation for r^ 4 ^ A (g), 



r (4)A (g)=U(g) + U(g)L A (g)f( 4 ) A (g) 



(86) 



where 



L A (g) = ^f(p)L A (p;g)f(p) 

p 



( Li(q) \ 

L A (q) Lf{q) L' A (q) 

-L>; A (q) L A (q) -If (q) 

V L' A (q) Lf(q) L A (q) J 



(87) 



The matrix elements L^- and L A with j 



1,2,3 vanish for symmetric form factors. The 
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other matrix elements are given by 

l k M = E [° A (p - i)° A (p + f) - pA (p - d fA o° + f)l fc (p) > 

p 

£ m(9) = E ^ " §) GA (P + 1) + ^ - f ^ + D] > 
p 

ift?) = -E{ Re [ GA (p+i) GA (-p+i)]- FA (p-i) FA (p+i)}^ 2 W' 

p 

~ L » = -E{ Re [ GA (p+i) GA (-p+i)]+ FA (p-f) FA (p+f)}/p(p)' 
L p a (i) = - E Im [ GA (p + D GA (-f + D] /» > 

p 

L' A (g) = 2^Re [G\p - 2 )F A (p + 2)] / c (p)/ p (p) , 
p 

£f (?) = 2 E Im - D fA (p + 1)] Up) Up) ■ (88) 



The system of linear equations Eq. (86) can be solved explicitly. The magnetic coupling 
function is decoupled from the others, so that M^\q) + M A (q) = {[M( )(g)] _1 - L A (g)}" 1 . 
Solving for the other coupling functions amounts to solving a linear 3x3 system. We do 
not write the explicit expressions here, but discuss only the singularity structure of the 
coupling functions. Singularities arise because the determinant d A (q) = det[l — U(g)L A (g)] 
vanishes at q = for A —> 0, if A A (k) remains finite, that is, in case of spontaneous 
symmetry breaking. For q — 0, the explicit solution for the coupling functions and their 
behavior for A — > was discussed in detail in Ref. [2H For small q ^ 0, one can expand 
d(q) = d A=0 (q) = do + d\q^ + c?2q 2 + . . . , where do oc Ao for small Ao, while d\ and d 2 remain 
finite for Ao — > 0. Expanding all coefficients to leading order in go and q, one obtains the 
following singular coupling functions: 

1 



$(g) oc 
P"(q) oc 



do + diql + d 2 q 2 ' 

go 

do + di9o + ^q 2 



X"(q) oc " u - - (89) 

w d + d x ql + d 2 q 2 V ' 

for A = 0. The other coupling functions, C (q), M A (q), A A (q), and the real part of X A (q) 

remain finite for A — > 0, Ao — > 0, q — > 0. 

The divergence of the vertex in the phase fluctuation channel represented by the coupling 

function $ A (g) reflects the Goldstone mode associated with the spontaneous breaking of the 
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U(l) symmetry. The Goldstone theorem, which guarantees the existence of this mode, is 
obviously respected by the RPA. A less familiar interesting result of the above calculation 
is the divergence of the (3+l)-interaction represented by the coupling function X A (q). At 
q = this interaction describes pair annihilation (or creation) combined with a forward 
scattering process. 



VI. ATTRACTIVE HUBBARD MODEL 



In this section we compute the flow of the Nambu vertex and the gap function for the 
two-dimensional attractive Hubbard model as a prototype of a spin-singlet superfluid. 
The Hubbard model describes interacting spin-^ lattice fermions with the Hamiltonian 

H = J2 *u4cj<r + U Wjtnu > (90) 

ij j 

where C; CT and c- xa are creation and annihilation operators for fermions with spin orientation a 
on a lattice site i. For the attractive Hubbard model the interaction parameter U is negative. 
The hopping matrix is usually short-ranged. We consider the case of nearest and next- 
to-nearest neighbor hopping on a square lattice, with amplitudes —t and — t r , respectively, 
yielding a dispersion relation of the form 

e(k) = —2t (cos k x + cos k y ) — At' cos k x cos k y . (91) 

The ground state of the attractive Hubbard model is a spin-singlet s-wave superfluid 
at any filling factor.^ For if — the superfluid order is degenerate with a charge den- 
sity wave at half-filling (only). The attractive Hubbard model has been studied already 
in several works both at zero and finite temperature by resummed perturbation theory 
(mostly T-matrix) GSEH quantum Monte Carlo (QMC) methods,'^®' and dynamical mean- 
field theory™! 



A. Regularization and counterterm 

The renormalization group flow is governed by the scale dependence of the regularized 
bare propagator, which we choose to be of the following form 



_! / ^ -e(k)-^ A (k) + i? A (fc ) A 



K(k)]- L = I ^ u SViV " s ViV ' " V " VJ "° I , (92) 

~ ' A ik + Z(k) + ^ A M + R A (k Q ) 1 
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with £(k) = e(k) — //. The regulator function 

R A (k ) = tsgn(k )^/k 2 + A 2 - tk (93) 

replaces frequencies ko with \ko\ <C A by sgn(A;o)A and thus regularizes the Fermi surface 
singularity of the bare fermionic propagator. The (real) bare gap A induces symmetry 
breaking and regularizes the Goldstone mode singularity forming in the effective interaction 
below the critical scale A c . Instead of linking the flow of Ao to the fermionic cutoff scale A 
by defining a A-dependent A A , we found it more convenient to keep Ao fixed until A has 
decreased to 0, and send Ao to zero afterwards. The equations for the latter flow are obtained 
simply by replacing A-derivatives by derivatives with respect to Ao- The counterterm <5£ A (k) 
is linked to the normal component of the self-energy by the condition 

^[^ A (k) + E A (0,k)] = 0, (94) 

such that the Fermi surface remains fixed during the flow. Since there is a contribution 



proportional to d\5£ to the scale derivative of E , solving Eq. (94) for d\5£ amounts to 
solving a linear integral equation. 



B. Parametrization 



We now specify the approximate parametrization of the self-energy and the interaction 
vertex. Due to the local bare interaction and the pairing instability occurring in the s- 
wave channel, the momentum dependence of the normal and anomalous self-energy can 
be expected to be weak, at least at weak coupling. Perturbation theory^ and previous 
functional RG calculations^ showed that this is indeed the case. We therefore discard the 
momentum dependence of the self-energy, keeping however the frequency dependence. The 
latter is treated numerically by discretizing E A (/co) and A A (ko) on a suitable grid. The 
counterterm <5£ A is then also momentum independent and can be interpreted as a shift of 
the chemical potential. 

The interaction vertex is fully described by the coupling functions C A fc , (q) etc. introduced 
in Sec. IV, where singular momentum and frequency dependences have been isolated in one 
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variable q. We now approximate these functions by the following ansatz: 

C A k ,{q) = C\q) g A (k ) g A (k' ) , 
M A k ,(q) = M A (q) g A (k ) g A (k' ) , 

a£M = A A (q) g A (k ) gt(K) > 
= * A (?) gi(h) , 
p®(q) = p"Hq) gi(ko) g^K) , 

X&b) = X' A (q)g A (k )g A (k' ), 

X&b) = X"\q) g A (k ) g A (k' ) . (95) 

The vertex thus assumes the form of a collection of boson-mediated interactions with bosonic 
propagators coupled to the fermions via fermion-boson vertices g A . The latter are normalized 
to one at zero frequency (k = 0). The momentum dependence on k and k' has thus 
been neglected, and the dependence on k Q and k' has been factorized. For the attractive 
Hubbard model, dependences on k and k' are generated only at order U 3 , and can thus 
be expected to be weak at least at weak coupling. Neglecting the dependence on k and k' 
implies a restriction to s-wave symmetry in charge, magnetic and pairing channels. As a 

A S .A. A S .A. 

consequence, all triplet components vanish, such that A£ k ,(q) = A k ' k , (q), <& kk ,(q) = § kk , (q), 



and X kk ,(q) = X kk A (q), and in the matrix V PP,A , Eq. (66), only four elements are non-zero. 
Compared to an exact decomposition of the coupling functions as in Eqs. (69) and (71), the 
sum over basis functions is replaced by just one term in the above ansatz. Due to time- 
reversal and exchange symmetries there is no contribution to W k {j(q) of that form. We have 
allowed for four distinct fermion-boson vertices g A , g^, g A i an d g A - The factorization of the 



coupling functions is similar to the factorization Eq. ( 79 ) obtained for separable interactions 
in RPA. Instead of parametrizing the fermion-boson vertices in the pairing channel by a 
single function g A , we now distinguish between g A and g A . It turns out that they differ only 
slightly. The imaginary part of the pairing coupling function P kk i(q) has little impact on the 
flow. Instead of introducing another fermion-boson vertex for that quantity, we approximate 
its dependence on k and k' by g A . The frequency dependence of the fermion-boson vertices 
g A {ko) is treated numerically by discretization. 

The parametrization of the "boson propagators" C A (q) etc. requires special care, to cap- 
ture the singularities. We first consider the amplitude and phase channel. For small q 
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the functions A A (q) and $ A (g) behave as [A A (q)]~ l = -m A - Z A q 2 - Z A q 2 + ... and 
[$ A (g)]- 1 = - TO A - Z A q 2 - Z^ql + . . . , where m A ->• for A < A c , A -»■ 0. Actually the 
regulator function can also generate contributions of order \q \, which disappear again as 
A — > 0. To deal with this technical complication, and to achieve an accurate parametriza- 
tion also at larger values of go and q, we parametrize A A and <3? A by two scale-dependent 
functions, 

[A^l)]- 1 = -mt(q ) - e A (q) , 

[* A (?)] _1 = -^(%)-e A (q), (96) 

where e A (0) = e A (0) = 0. The (discretized) momentum and frequency dependences of these 
functions are determined from the flow. The above ansatz with functions of one (go) and 
two (q = (q x ,Qy)) variables reduces the numerical effort compared to a discretization of an 
arbitrary function of g = (qo,q x ,Qy)- Tests within RPA indicate that it describes the full 
functions sufficiently well. In particular, the behavior at small go and small q is captured 
correctly. The imaginary parts P" A (q) and X" A (q) are odd functions of go- This and the 



expected singularity structure [see Eq. (89)] motivate the ansatz 



P"\q) =r q ° 



m p»(lo) + e A „(q) ' 
- V "" V/i = ~ a i \~T a / \ ' (97) 



The parametrization of C A (g), M A (g) and X /A (g) is slightly more complicated, because 
at small q these functions cannot be represented as a sum of a frequency and a momentum 
dependent function. We therefore distinguish the cases |q| < g max and |q| > g max with a 



suitably chosen g max - For |q| > g max we make an additive ansatz analogous to Eq. (96), 

[C A (g)] _1 = -m A (g ) - e A (q) , 
[M A (g)r X = -mi (g ) - e A (q) , 

[X'^q)]- 1 = -m A (g )-e A (q). (98) 

For small q, the q-dependence is increasingly isotropic, except for the special case where 
the Fermi surface touches van Hove points (which we exclude). Hence, for |q| < g max we 
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approximate the momentum dependence as isotropic, 

C\q) = C A (g ,|q|), 
M A (g) = M A (g ,|q|), 

X'\q) = X' A (g ,|q|), (99) 

reducing the number of variables again to two. To avoid a discontinuity at g max , we connect 
the two regimes in momentum space by a smooth partition of unity instead of step functions. 



C. Flow equations 

The flow equations for the scale-dependent functions parametrizing the self-energy and 
interaction vertex are obtained by projecting the flow equations for the self-energy and the 
coupling functions on the simplified ansatz. Dependences on the fermionic momenta k, k' 
generated by the flow are eliminated by a Fermi surface average (but not the q-dependence, 
of course). This corresponds to keeping only the leading (in power-counting) term in an 
expansion around the Fermi surface, and averaging over the momentum dependence along 
the Fermi surface, which is in line with the pure s-wave ansatz for the interactions. 

For the self-energy, we project on the momentum-independent ansatz by averaging the 
flow Eq. Q over the Fermi surface as follows: 

^E A (A; ) = (rhs A (A;o,k)) keFS , (100) 

where rhs A (/c ,k) stands for the right hand side of the flow equation (in Nambu matrix 
form), and (. . .)keFS denotes a Fermi surface average. Momentum dependences of the self- 
energy perpendicular to the Fermi surface are marginal in power-counting^ and lead to a 
(finite) renormalization of the Fermi velocity. However, they are quantitatively small in the 
attractive Hubbard model, at least at weak coupling and away from van Hove points, and 
have thus little influence on our results. 

The flow equations for the coupling functions C£ k ,(q), . . . , X£ k ,(q) were derived in Sec. IV. 
Inserting the ansatz for the interaction vertex on the right hand side of these equations 
yields several terms which can be represented by Feynman diagrams of the form plotted in 
Fig. 6. We recall that the point-like vertex represents the bare (here Hubbard) interaction, 
the wiggly line any coupling function contributing to V PH ' A , and the double line coupling 
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functions appearing in V PP ' A (only Mj} k ,(q) in the absence of triplet terms). Note that 
the terms in Fig. 7 are redundant since the complete set of coupling functions is already 
captured by V PH,A . The Nambu index sums on the right hand side of the flow equation for 
the coupling functions can be transformed to a more convenient form by representing the 
vertex and the propagator product in the Pauli matrix basis defined in Appendix A and 
used already in Sec. V. 

Some of the contributions, having the form of vertex corrections or box diagrams, generate 
dependences on k and k' which are not allowed for in our ansatz. These dependences are 
projected out by a symmetrized Fermi surface average. We discuss the procedure for the 
charge coupling function as a prototypical example, which can be extended directly to all 
other cases. The flow of the projected charge coupling function is given by 

d [C\q) g^ko) = (rhs A (M';g)> k±f , k , ±feFS 

= ihs A {k ,k' ;q) , (101) 



dA 



where rhs A (/c, k'; q) denotes the complete right hand side of the flow equation for C£ k ,, 



Eq. (|40j), and 

(• • -)k±|,k'±f 6FS = | X] (• • -)k+4,k'+e'f 6FS (102) 
e,e'=±l 

is a symmetrized Fermi surface average. The latter averages the four possible ways of 
integrating k and k' under the constraint that two of the four external momenta k ± | 
and k' ± | of the vertex lie on the Fermi surface. For q = 0, corresponding to the forward 
scattering and Cooper channels, this becomes a Fermi surface average with all four momenta 
on the Fermi surface. For q ^ 0, the set of momenta k satisfying both k + 3 e FS and 
k — | £ FS is limited to few points, except for special nesting vectors in case of a nested 
Fermi surface. The Fermi surface average picks up the s-wave component of the dominant 
processes near the Fermi surface. Indeed, the momentum dependence perpendicular to the 
Fermi surface is irrelevant in power-counting.^ Note, however, that we do not discard the 
dependence on q. That dependence becomes important due to the formation of singularities, 
which invalidate the weak-coupling power-counting. 

The projection on the form factors in the channel decomposition could also be carried out 
by integration over the entire Brillouin zoneP^H However, for our simple ansatz with only 
one (momentum-independent) form factor it is better to approximate the vertex by its Fermi 
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surface average instead of a Brillouin zone average, to capture the dominant contributions 
at low energy scales. We checked this in some test cases by explicit comparison of different 
projection procedures. 



The flow equation for the bosonic propagator can be extracted from Eq. ( 101 ) by setting 
ko = k' Q = 0. Since <? A (0) = 1 is independent of A, one obtains 

^ K C A (q)=Ths\0,0;q). (103) 

The functions m A (q ) and e A (q) parametrizing C A (q) for |q| > g max are extracted by evalu- 
ating [C A (q)]^ 1 at a fixed momentum q* as a function of qo, and at fixed frequency go = as 
a function of q, respectively. For q* we choose a momentum where C A (q) is peaked, where 
it yields the largest contribution. In the charge and magnetic channel this happens typically 
at finite momenta connecting antipodal Fermi points (2fcp-type momenta). 



The product rule for differentiation applied to the left hand side of Eq. (101) at k' = 
yields the flow equation for the fermion-boson vertex, 

d -9c(ko) = 7^ te A (fc , 0; q) - g A (k ) rhs>, 0; qj\ ^ . (104) 



dA *c v u, 



9=(0,q*) 



The flow equations for the other coupling functions M£ k ,(q) etc. are projected on the ansatz 
in the same way. The flow of the fermion-boson vertices in the pairing channel g A {ko) and 



g A (k ) is determined as in Eq. (104), with q* = 0. 

The initial conditions for the flow at Ao = oo are as follows. For the self-energy, coun- 
terterm and gap function the flow starts at S A ° = 0, 5£ A ° = 0, and A A ° = A . The coupling 
functions are initially zero, and the fermion-boson vertices are equal to one. Note that the 
coupling functions do not include the bare interaction. In a numerical evaluation, the flow 
starts at a large finite Ao- The self-energy receives a tadpole contribution of order one in 
the flow from Ao = oo to an arbitrarily large finite Ao, yielding S Ao — U/2 with corrections 
of order Aq 1 , and correspondingly <5£ A ° = —U/2. The error of order Aq 1 made by starting 
the flow at a (large) finite cutoff can be significantly reduced by using perturbative results 
at A as initial conditions instead of the initial values at A = oo. The coupling functions 



are then non-zero from the beginning such that Eq. (104) is well defined at Ao- 

We conclude this section with a few words on numerical aspects. More details can be 
found in Ref . US Momentum and frequency dependences were discretized on non-equidistant 
grids such that the resolution is higher at smaller frequencies and momenta. The positive 
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frequency axis and radial momentum dependences were discretized by around 30 points, 
and angular momentum dependences by 6 angles per quadrant in the Brillouin zone. The 
functional flow equations were thus replaced by a system of around 2000 non-linear ordinary 
differential equations with three-dimensional loop integrals on the right hand sides. The 
integrals were performed with an adaptive integration algorithm and the integration of the 
flow was performed with a third-order Runge-Kutta routine. Depending on parameters, the 
computation of a flow required between a day and a week on 20 CPU cores. 



D. Results 



We now present results for the effective interactions, the normal self-energy, and the gap 
function as obtained from a numerical solution of the flow equations. Most of the numerical 
results are obtained for a small fixed external pairing field A chosen two to three orders of 
magnitude below the mean-field gap Amf, but we also discuss some flows where Ao scales 
toward zero after the fermionic cutoff has reached A = 0. The Ward identity following from 
global charge conservation is enforced at zero frequency by projecting the flow, if not stated 



otherwise (for details, see Sec. VI D 3). Bare interaction strengths are chosen in the weak to 
moderate coupling range \U\/t = 1 — 4. In the following we set the nearest-neighbor hopping 
amplitude t — 1, that is, all quantities with dimension of energy are in units of t. 



1. Effective interactions 

We begin with results for the coupling functions, which describe the various effective 
interaction channels contributing to the the Nambu vertex. With our sign conventions 
negative coupling functions correspond to attractive effective interactions in the respective 
channel. 

The flow of effective interactions in the pairing channel is qualitatively similar to the one 
in RPA (see Sec. V). However, the critical scale and the size of the coupling functions is 
reduced by fluctuations. Typical flows for the amplitude and phase couplings at q = are 
shown in Fig. 8 for various choices of the external pairing field Ao- For U = —2 stable flows 
without artificial singularities could be performed for external pairing fields as small as three 
orders of magnitude below the size of the mean-field gap Amf, with a final phase coupling 
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FIG. 8: (Color online) Scale dependence of the amplitude and phase couplings A A and <3? A at q = 
for various choices of the external pairing field Ao- The Hubbard model parameters are t' = —0.1, 
U = — 2, n = 0.5 (quarter- filling) . 

$ A=0 (0) proportional to Aq 1 within numerical accuracy, as dictated by the Ward identity. 
The amplitude coupling A A (0) has a peak around A c , whose size increases strongly upon 
reducing Ao, while it reaches a finite value with a negligible dependence on the external 
pairing field at the end of the flow. 

The momentum and frequency dependence of A A (q) and $ A (g) around q = is shown 
for various choices of A in Fig. 9. For small momenta, the coupling functions are isotropic 
functions of q with a momentum dependence proportional to q 2 , for both finite A and 
A = 0. The frequency dependence is linear for small go at A > 0, but essentially quadratic 
for A = 0. The linear behavior at A > is caused by the frequency dependent regulator, 



Eq. (93), and thus disappears once the regulator has scaled to zero. The amplitude coupling 
A A=0 (q , 0) exhibits a tiny dip at q = 0. Overall, the qualitative momentum and frequency 
dependences of the coupling functions in the pairing channel do not deviate significantly 
from the behavior in RPA. This is also true for the imaginary part of the pairing coupling 
P" K {q) and the imaginary part of the anomalous (3+l)-coupling X" A (q), whose singular 
behavior at small momenta and frequencies is well described by 

-//A=0/„\ „ ry/A=0/_A Qo 



X" A=0 (q) oc P //A=u (g) cx 



(105) 



A + aql + bq 2 ' 
where a, b are positive constants. 

The charge coupling function C A (q) is generally negative at all stages of the flow. It 
thus renormalizes the bare attraction in the charge channel given by U to an enhanced total 
attractive interaction U + 2C A (q). This effect is captured already in RPA. The enhancement 
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FIG. 9: (Color online) Momentum dependence along the g x -axis (left) and frequency dependence 
(right) of the amplitude and phase couplings A A (q) and & A (q) for small momenta and frequencies 
at various stages of the flow. The Hubbard model parameters are the same as in Fig. 8 and 
A = 10- 3 . 

is usually small. However, for densities near half-filling and small values of if it becomes 
large at q = (0, Q) with Q = (71", 7r). For if = and half-filling, U + 2C A (0, Q) is degenerate 
with the pairing interaction U + P (0,0), reflecting the degeneracy of superfluidity with 
charge density wave order due to a particle-hole symmetry in this special caseP^ In Fig. 10 
we show the momentum dependence of the charge coupling function in the static limit go — 
at the end of the flow (A = 0) for two distinct choices of t' and n. The function exhibits 
pronounced peaks at incommensurate momenta situated at the Brillouin zone boundary. 
These peaks are present already in the bare polarization function (particle-hole bubble).^ 
They move toward (it, ir) and increase upon approaching half-filling for t' = 0. As a function 
of frequency, the size of C A (g) decays monotonically upon increasing |go|. 

The magnetic coupling function M A (g) receives contributions beyond RPA which change 
its behavior qualitatively. In Fig. 11 we show its momentum dependence in the static limit 
go = at the end of the flow for the same choices of t' and n as in Fig. 10. The coupling 
function is negative in most of the Brillouin zone, but it develops a pronounced positive 
peak for small momenta q. This peak is a pure fluctuation effect. In RPA, M A=0 (0,0) 
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FIG. 10: Momentum dependence of the static charge coupling function C °(0, q) at the end of 
the flow, for t' = —0.1, n = 0.5 (left) and t' = 0, n = 0.78 (right). The Hubbard interaction is 
U = —2 in both cases. 




FIG. 11: Momentum dependence of the static magnetic coupling function M A=0 (0, q) at the end 
of the flow, for t' = —0.1, n = 0.5 (left) and t' = 0, n = 0.78 (right). The Hubbard interaction is 
U = —2 in both cases. 

vanishes due to the pairing gap. Since the amplitude of the coupling function is small, the 
total interaction in the magnetic channel — U + 2M A=0 (q) is dominated by the bare Hubbard 
interaction and remains positive for all momenta. In Fig. 12 the frequency dependence of 
M A (q Q , 0) is shown at various stages of the flow. The positive peak at go = develops at 
and below the critical scale A c and is foreshadowed by a finite frequency peak for A near A c . 
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FIG. 12: (Color online) Frequency dependence of the magnetic coupling function M A (go,Q) at 
various stages of the flow above (left) and below (right) the critical scale for pairing. The model 
parameters are t' = —0.1, U = —2, and n = 0.5. 

The flow is non-monotonic and M A (go, Q) exhibits a sign change for small go, but eventually 
M A=0 (go, Q) is positive for all frequencies. A similar sign change at finite go and a pronounced 
finite frequency peak has been observed previously in the charge coupling function for the 
repulsive Hubbard model in the symmetric regime (A > A e )J 45 * 48 l 

The real part of the anomalous (3+l)-coupling function X /A (g) is relatively small. Its 
singularity at the critical scale A c is considerably broadened by fluctuations (beyond RPA). 
Nevertheless, its influence on the flow of the self-energy and the other coupling functions is 
important. Neglecting the (3+l)-coupling would lead to artifacts like non-monotonic flows of 
$ A (0) even for small interactions U. While the imaginary part of A A (g) depends strongly on 
A for small q and A, the real part does not. In Fig. 13 we plot the momentum dependence 
of the static (3+l)-coupling function at the end of the flow for the same choices of t' and n 
as in Figs. 10 and 11. Note that the imaginary part of the static (3+l)-coupling function 
vanishes, that is, X A (0,q) is real. 

We now turn to the fermion-boson vertices, whose frequency dependence is plotted in 
Fig. 14. Note that the vertices are even functions of go which are normalized to one at 
g = by definition. The frequency dependence of the vertices is quite weak. However, 
the frequency dependence of the vertices in the pairing channel contributes significantly 
to the frequency dependence of the gap function A A (fco) and also to the flow of $ A . The 
normal self-energy and the other coupling functions are only weakly affected by the frequency 
dependence of the fermion-boson vertices. The magnetic vertex exhibits a small peak at low 
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FIG. 13: Momentum dependence of the static (3+l)-coupling function X A=0 (0,q) at the end of 
the flow, for t' = —0.1, n = 0.5 (left) and t' = 0, n = 0.78 (right). The Hubbard interaction is 
U = —2 in both cases. 
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FIG. 14: (Color online) Frequency dependence of the fermion-boson vertices at the end of the flow 
(A = 0). Left: amplitude and phase vertices. Right: charge and magnetic vertices. The model 
parameters are t' = —0.1, U = —2, and n = 0.5. 

frequencies which develops at scales A < A c and is therefore related to pairing fluctuations. 
2. Normal self-energy and gap function 

At weak to moderate interactions the ground state of the attractive Hubbard model is 
superfluid with Cooper pairs made of weakly renormalized quasi particles. Quasi particle 
renormalization occurs already at scales above the pairing scale A c and is described by the 
normal self-energy. The momentum dependence of the self-energy is weak for the choice of 
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FIG. 15: (Color online) Frequency dependence of the imaginary part of the normal self-energy for 
t' = —0.1, n = 0.5 (quarter-filling) and various choices of U at the end of the flow. The result 
labeled as "no g(ko)" is obtained with constant (frequency-independent) fermion-boson vertices. 

parameters considered in this work and we will present only results for the Fermi surface 
average E A (&o) = (£ A (/co, k))k e ps • In Fig- 15 we show results for the imaginary part of 
E A (/co) as a function of frequency at the end of the flow (A = 0). We plot only the positive 
frequency axis since Im£ A (— k ) = — ImS A (/c )- The real part (not plotted) of S A (fc ) is 
an even function of k with a negative peak at k Q = that decays monotonically to the 
Hartree term Un/2 with increasing \ko\. The overall shape of the self-energy is the same for 
all interaction strengths, only the size increases with \U\. 

The slope of Im£ A (fc ) at k = yields the quasi particle weight Zf as 



Z f = Zf=° ranges from Z f = 0.96 for U = -1.5 to Z f = 0.87 for U = -3. Although the 
normal self-energy is fairly small and the quasi particle weight is only slightly suppressed 
for small to moderate interactions, it has nevertheless a significant impact on the size of the 
pairing gap. 

Flows of the gap A (/co) at ko = are shown in Fig. 16 for various choices of U. The 
small external pairing field Ao increases to much larger gaps at scales near and below the 
critical scale A c , where v4 A (0) has a peak. The edge of the gap flow at A c becomes sharper for 
smaller A . The gap at the end of the flow assumes values close to A c . The scale dependence 
for A < A c obeys approximately 




(106) 



A A (0) w ^A 2 C - A 2 , 



(107) 
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FIG. 16: (Color online) Scale dependence of the gap A A (&o) at fco = for t' = —0.1, n = 0.5 and 
various choices of U and Ao- For U = —2 and —3, the mean-field scale-dependence \/A;? — A 2 , 
with A c determined from the peak of t4 a (0), is shown for comparison. 

with increasing accuracy for smaller values of U and Ao- In mean- field theory this relation 
is exact for A — > 0, as one can easily see by writing down the gap equation in the presence 



of the infrared regulator Eq. (93). For U = —2 the gap flow lies almost on top of the square 
root function Eq. (107) for small A , while for stronger attractions deviations become visible. 
In particular, the final gap A A=0 becomes clearly larger than A c . 

The flow in Fig. 16 was obtained with frequency dependent effective boson propaga- 
tors and fermion-boson vertices, and a frequency dependent normal self-energy and gap as 



described in Sec. |VI B| Comparing with results obtained by discarding the frequency depen- 
dence of some of these quantities, one finds that only the frequency dependence of the boson 
propagators and of the imaginary part of the normal self-energy have a substantial impact 
on the size of A A (0). The feedback of the other frequency dependences on the gap at k = 
is small. 

The critical scale and the final gap are strongly reduced compared to their mean-field 
values A^ F and Amf, respectively, mostly due to fluctuations above the critical scale. In 
Fig. 17 we plot the ratio A/Amf with A = A A=0 (0) as a function of U for t' = —0.1 and 
n = 0.5. The lower curve was obtained by a simplified static parametrization of the vertex 
and self-energy, where all frequency dependences where neglected. Notably the reduction 
increases at weaker interactions and does not extrapolate to one for U — > 0. This is actually 
the expected behavior. In the weak coupling limit the gap A has the same exponential 
[/-dependence A oc e _b// ' C7 ' with a (density-dependent) constant b as in mean-field theory. 
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FIG. 17: (Color online) Gap ratio A/Amf as a function of U as obtained from the flow with 
frequency dependent (upper curve) and static (lower curve) vertices and self-energies. 

However, the prefactor of the BCS mean-field formula is reduced by fluctuations, as first 
noted for the transition temperature in three-dimensional superconductors by Gorkov and 
Melik-BarkhudarovpSl The reduction factor in the weak coupling limit can be computed 
by second order perturbation theory. ^ * 50 * 51 ! For the parameters used in Fig. 17 one finds 
A/Amf — > 0.3 for U — » 0. 46 Both curves in Fig. 17 should tend to that value, since the 
flow captures the perturbative contributions. However, we cannot reach the limit U — > 
numerically. It is hard to compute the gap from a numerical solution of the flow equations 
for smaller interaction strengths than those shown, since A c and A decrease exponentially. 

For strong attractions U the attractive Hubbard model can be mapped to a Heisenberg 
model in a uniform magnetic field The gap ratio A/Amf thereby translates to the ratio 
between the staggered magnetization m s and the corresponding classical result mf. From 
numerical results for that raticJ^ one can infer that the gap ratio in the strongly attractive 
Hubbard model is 0.6 at half-filling and even larger away from half-filling. The observed 
increase of A/Amf with increasing \U\ is therefore consistent with the expected trend. 
Similar values for A/Amf but with a less pronounced [/-dependence have been obtained in 
an earlier fRG study with a simpler parametrization of the vertex.^ 

We now discuss the frequency dependence of the gap function. In Fig. 18, A A=0 (k ) is 
plotted as a function of frequency for U = —2, i! = —0.1, and n = 1/2. Results obtained by 
computing the gap from a projected flow obeying the Ward identity at k = are compared 
to results where the frequency dependence of the gap is computed directly from the Ward 
identity (A\yi), contrasting also calculations with and without frequency dependent fermion- 
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FIG. 18: (Color online) Frequency dependence of the gap at the end of the flow for various ap- 
proximation schemes. The inset shows the gap at small frequencies fco < 1. The model parameters 
are U = -2, t' = -0.1, and n = 1/2. 

boson vertices g(ko). Note that the results discussed so far were all obtained by enforcing 
the Ward identity only at k = 0. Unlike the frequency dependence of the normal self- 
energy, the frequency dependence of the gap is strongly affected by the frequency dependent 
renormalization of the fermion-boson vertices. Neglecting it leads to a very weak (A) or 
almost no (Awi) frequency dependence. The gap A(k ) computed by projecting the flow on 
the Ward identity at ko = exhibits a shallow finite-frequency minimum, which is probably 
an artifact of the approximations, associated with a (slight) violation of the Ward identity 
at finite frequencies. Awi(fco) nas a minimum at k = 0. A qualitatively similar frequency 
dependence of the gap is also captured by the T-matrix approximation.^ 

3. Ward identity 

For real gaps A and A A the Ward identity, Eq. rt8J), can be simplified to 



A A 



(fc)-Ao(AO = -5> (A/) \G\k>)G\-k>) + {F\k')f 



k' 




(108) 
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Expressing V A and W A by the coupling functions introduced in the channel decomposition 
(Sec. IV), the identity can be written as 



&\k) = -J2Mk') G\k')G\-k') + (F\k')) 2 $a,( ) + O(A 



(109) 



for A — > 0. The first term on the right hand side is of order one, since 3^/(0) oc A 1 
for small A . With the approximate parametrization for the Hubbard model described in 



Sec. VI B, the combination of interaction terms on the right hand side of Eq. (108) can be 



written as 



V A (k, —k, —k', k') - W A (k, k', k', k) 
+ l -A\k'~k) [g A (p, 
+ C\k'-k) [g A {p, 



U + $ A (0)g A (k o )g A (k' o ) 
\$ A (k' - k) [g A (p )} 2 
3M A (k> - k) [gi( Pl 



(110) 



where p$ = (ko + k' )/2. For a small constant Ao and a momentum-independent gap function 
A A (/c ), the Ward identity then assumes the form 

A A (M = -J2 A ° [G A (k')G A (-k') + {F A (k')) 2 ] ®\0)g A (k )g A (k' ) + O(A ) . (Ill) 

k> 

The most important consequence of the Ward identity is the divergence of the phase 
coupling $ A (0) in the limit Ao — > for A < A c , reflecting the massless Goldstone boson 
associated with spontaneous symmetry breaking. The truncated flow equations do not 
obey the Ward identity exactly, and $ A (0) deviates from the expected behavior oc Aq 1 
for small A . For small U the deviations are tiny. For example, for t' = —0.1, n = 1/2, 
and U = —2, the product A $ A=0 (0) is almost constant down to fairly small values of A , 
before it increases and finally diverges at a finite Ao of the order 10~ 5 , which is four orders of 
magnitude smaller than A A=0 . The same behavior was observed already in more pronounced 
form in Ref . [23l 

The violation of the Ward identity can be quantified by comparing the gap A^ G computed 
from its flow equation to the gap A^j required by the Ward identity. The latter is computed 



from Eq. (108) by inserting the coupling functions as determined from the flow on the right 
hand side. In Fig. 19 we plot the difference A^ G — A^j, divided by A A ^ G U /L , as a function 
of the scale in units of A c . One can see that the violation builds up gradually at scales 
around A c . The normalized difference (A^ G — A A vl )/A A lG increases rapidly from U = —2 
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FIG. 19: (Color online) Violation of the Ward identity as a function of the scale A for t' = —0.1, 
n = 1/2 and various values of U and Ao- The gap A^ G determined from the flow equation is 
compared to the gap A^j determined from the Ward identity. 

to U — —3. For A > A c it is roughly proportional to U A . For A < A c a pronounced 
Ao-dependence appears. For much smaller values of Ao than those shown, A^ G — A^i can 
turn negative, which is related to an artificial divergence of $ A at a finite A . We observed 
similar [/-dependences also for other hopping parameters and densities. On general grounds 
one would expect a violation of the Ward identity of order U 3 at weak coupling, even if 
the one-loop flow was carried out without additional approximations.^ 3 ^ 1 The above results 
suggest that the violation sets in only at order U A , or contributions of order U 3 have very 
small prefactors. 

In Fig. 20, ( A^ G — A^ : ) / A^_ G is plotted as a function of A for a fixed set of parameters, to 
compare the performance of different approximations. The graph labeled "dynamical" was 
obtained by using the frequency dependent parametrization of the vertex and the self-energy 



as described in Sec. VI B The graph "no g(fco)" was computed with constant (unity) fermion- 



boson vertices, and the graph "static" by discarding all frequency dependences. The lowest 
curve labeled "coord, proj." was computed with the dynamical parametrization and the 
Ward identity enforced by "coordinate projection" (as described below). The latter obeys 
the Ward identity by construction, up to small discretization errors. Taking the frequency 
dependences into account obviously reduces the violation of the Ward identity significantly. 

Even for the most accurate parametrization of the vertex, the Ward identity is not fulfilled 
by the truncated flow, as generally expected. 1 27 1 46 1 The deviations are small for weak interac- 
tions but increase rapidly with \U\. Violating the Ward identity spoils the singular infrared 
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FIG. 20: (Color online) Violation of the Ward identity as a function of A for U = —2, t? = —0.1, 
n = 1/2, and Ao = 10 -3 . Results from the flow equations with static and two distinct (with and 
without frequency dependent fermion-boson vertices) dynamical parametrizations of the vertex are 
compared to the result from the Ward identity projected flow. 



behavior of the coupling functions associated with the massless Goldstone boson for A — > 0. 
Even worse, it leads to artificial singularities which prevent one from carrying out the flow 
down to A — > and Ao — > 0. In the results presented in the preceding sections we have 
therefore enforced the Ward identity by using a coordinate projection procedure, devised 
for the numerical solution of systems of ordinary differential equations with constraints.^ 
The flowing quantities are thereby projected on the manifold spanned by the constraint 
(Ward identity) in a way that the projected solution stays as close to the solution of the 
flow equations as possible, while deviations from the constraint are damped exponentially. 
In practice, we have enforced the Ward identity only at zero frequency (ko = 0), to reduce 
the numerical effort.®' This has little effect on absolute values of results, but leads to the 



slightly artificial frequency dependence of the gap at low frequencies discussed in Sec. VI D 2 



4- Aq- flow and singularities 

We finally take a closer look at the singularities of the vertex in the limit Ao — > 0. 
In particular, we complement the numerical results for the Hubbard model by qualitative 
analytical estimates which are generally valid for fully gapped singlet superfluids. 

To this end, we assume that the fermionic cutoff has already been removed (A — > 0), and 
we analyze the flow as a function of a decreasing pairing field Aq. In Fig. 21 we show the 
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FIG. 21: (Color online) A -flows of A(0), $(0) and A(0) for an initial value A = 0.005. Results 
obtained for fixed smaller values of Ao are shown for comparison (symbols). The model parameters 
are t' = -0.1, n = 0.5, U = -2. 

flow of A(0) = A A=0 (0), $(0) = $ A=0 (0) and A(0) = A A=0 (0) as a function of A , with an 
initial value A = 0.005. Results obtained for some fixed smaller values of A are shown 
for comparison. The numerical computation of the A-flow becomes increasingly difficult at 
smaller Ao- Furthermore, there are systematic deviations between the results obtained from 
the Ao-flow and those computed at fixed Ao- These may be related to divergencies in box 
diagrams for A — > 0, which can and must be treated by a flow starting at an initially finite 
Ao, as we discuss in the following. 

In a fully gapped superfluid, the fermionic propagator is regularized by the pairing gap. 
However, the interaction vertex develops a singularity associated with the emergence of a 
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Goldstone boson. In particular, the phase coupling function has the singular form 

*(g)«- A +z L + Zq2 . (H2) 

for small q = ((?Q,q), where Z qo and Z q are positive constants. This singularity is dictated 
by the Ward identity. Related singularities occur also for the imaginary parts of the pairing 
and anomalous (3+l)-coupling functions P"(q) and X"(q), respectively, but their impact is 
reduced by a numerator proportional to go- The divergence of <&(q) for q — > 0, A — > is 
integrable in (2+1) dimensions. Hence, self-energy and vertex correction (Fig. 4) contribu- 
tions involving integrals over <&(q) remain finite for Ao — » 0. However, box diagrams (Fig. 5) 
yield contributions involving an integral over products of two phase coupling functions, 

^Ao [Gs lS2 (p ~ q/2)G S3S4 (p + q/2)] $(p - k)$(k' - p) , (113) 

p 

where G is the propagator and $ the phase coupling function for A = 0. The Ao-derivative 
of the propagators is finite for Ao — > 0, but for k = k' the singularities of the phase coupling 

— 1/2 

functions coalesce and the integral diverges as A . It is therefore not possible to set A 
to zero before the fermionic cutoff has been removed. The A -flow is however well defined 
and integrable. 

Hence, the singularities associated with the Goldstone mode do not lead to divergencies 
in other channels. In this respect the one-loop flow analyzed in this work is qualitatively 
similar to the RPA. The fluctuation effects beyond RPA yield only finite renormalizations. 
On the other hand, it is known from the theory of interacting bosons that the phase mode 
does lead to a singular renormalization of the amplitude mode.^In a renormalization group 
theory of fermionic superfluids with auxiliary boson fields representing the order parame- 
ter fluctuations, this effect appears already at one-loop level.^ The singular contributions 
involve scale derivatives acting on the boson propagators. In the purely fermionic renormal- 
ization group (without auxiliary bosons), analogous singular contributions appear only at 
the two-loop levelpSl 



VII. CONCLUSION 



We have analyzed ground state properties of a spin-singlet superfluid including fluctua- 
tions on all scales via a fermionic functional renormalization group flow in a formulation that 
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allows for symmetry breaking. The flow equations were truncated in a one-loop approxima- 
tion with self-energy feedback. Spin rotation invariance and discrete symmetries were fully 
exploited to simplify the structure of the Nambu two-particle vertex. To parametrize the sin- 
gular momentum and frequency dependences of the effective interactions, the Nambu vertex 
was decomposed in charge, magnetic, and various normal and anomalous pairing channels, 
which are all mutually coupled in the flow. We have shown that the channel decomposed 
one-loop flow equations are equivalent to the RPA for the vertex and to mean-field theory for 
the gap function, if only direct Nambu particle- hole contributions are taken into account 
The crossed particle-hole and the particle-particle (in Nambu representation) contributions 
to the complete one-loop flow thus capture fluctuations beyond mean-field theory and RPA. 

We have evaluated the flow equations for the two-dimensional attractive Hubbard model 
as a prototype of an interacting Fermi system with a spin-singlet superfluid ground state. 
The dominance of s-wave terms in the effective interactions in that model allows for a rel- 
atively simple parametrization. The global U(l) Ward identity relating the vertex to the 
gap function is violated by the one-loop truncation. The deviations are very small for a 
weak attraction, but increase rapidly for stronger interactions. To maintain the singular- 
ity structure associated with the Goldstone boson, the flow was therefore projected on the 
Ward identity, analogously to evaluating a differential flow in the presence of a constraint. 
We have computed the effective interactions in the charge, magnetic, and pairing channels, 
including anomalous (3+l)-interactions describing pair annihilation (or creation) combined 
with a one-particle scattering process. Unprecedented comprehensive results on the mo- 
mentum and (imaginary) frequency dependences of the effective interactions were obtained 
and discussed. The singularities in the pairing channels generated by the one-loop flow are 
qualitatively similar to the RPA, and are to a large extent fixed by the Ward identity. The 
effective magnetic interaction develops a low-frequency small-momemtum peak which is a 
pure fluctuation effect. There are also significant quantitative fluctuation effects which are 
captured by the one-loop flow. In particular, the gap is strongly reduced compared to the 
mean-field value, with a stronger reduction at weaker interactions, as expected from pertur- 
bative and numerical results. The expected divergence of the superfluid amplitude mode in 
the low-energy limit is not captured by the one-loop truncation. This effect appears only at 
the two-loop level in the fermionic renormalization group flowP 

Besides the channel decomposition of the vertex for a system exhibiting spontaneous 
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symmetry breaking, there are two other noteworthy technical upshots of our work, which 
may be picked up in future calculations. First, we have found that an accurate discretization 
of both momentum and frequency dependences is computationally feasible^ and has several 
advantages compared to the usual strategy of an ansatz with a small number of scale- 
dependent coefficients. In particular, one avoids problems with momentum or frequency 
derivatives which are necessary to extract the flow of such coefficients. Second, we have 
shown that a symmetry breaking field can be used as a convenient flow parameter, which 
regularizes the flow at the critical scale and allows for a controlled treatment of infrared 
divergences associated with the Goldstone boson. 
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Appendix A: Pauli matrix basis 

It is often convenient to represent the Nambu vertex in a basis spanned by tensor products 
of Pauli matrices and the unit matrix. The Pauli matrices 

r (D ) r (2) > r (3) 

and the unit matrix 

r (o) f orm a basis in the vector space of complex 2x2 matrices. The tensor products r^'0r^ 
form a basis in the space of complex 4x4 matrices. The components of the Nambu vertex 
in this basis are obtained as 

f$ A (*!, k 2 , k 3 , h) = \Yl r iiW£ ^si 3 sM, fc, h, h) . (Al) 

Si 

The inverse basis transformation is given by 

rg^(fci, k 2 , k 3 , h) = 1 -J2 t&tW f$V, k 2 , k 3 , h) . (A2) 

3,3' 

The matrix formed by the components r^/ A is denoted as f ( 4 ) A . The tilde is used to 
distinguish this and other matrices represented in the Pauli basis from matrices in the 



Nambu index basis defined in Eq. (31 ) 



The flow equations for the coupling functions parametrizing the channel decomposed 
Nambu vertex can be derived most conveniently in the Pauli matrix basis. Since the complete 
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set of coupling functions is contained in the particle-hole contribution to the vertex, their 



flow is determined by the flow equation Eq. (60). Transformed to the Pauli matrix basis, 
the equation reads 



A-" 



where 



p hi' 



(A3) 



(A4) 



The decomposition Eq. (59) of the Nambu vertex can be written in the Pauli matrix basis 



with momentum variables fcj.,4 = k ± q/2 and A; 2) 3 = k' =)= q/2 as 

fg) A (*, k'; q) = U, f (k,k';q) + V^ A (k,k';q) 

t>PH',A ( k+k'-q k+k'+q . i, ,, 

v jj> y 2 ' 2 ' K K 

+ V^ A (^, k -^;k + k') . 



(A5) 



Note that V?, H ' A is defined by transforming V.^?'^ aj with the first two Nambu indices ex- 
changed to the Pauli matrix basis. The functions Lf-,(p; q) are given by products of normal 



and anomalous propagators, 



T A 

^0( 

r A 



T A 

-^03 



T A 
L A 
L 
L 



T A 



J 23 



(p;q 
(p;q 
(p;q 
(p;q 
(p;q 
(p;q 
(p;q 
(p;q 
(p;q 
(p;q 



R e [G A (p_)G A (p + )]+ F A (p_)F A (p + ) , 
zF A (p_)ImG A (p + ) + zImG A (p_)F A (p + ) = L A (p; q) , 
iReG A (p_)F A (p + ) -iF A (p_)ReG A (p + ) = -L A (p;q) 
iIm[G A (p^)G A (p + )]=L A {p;q), 
-Re{G A (p^)G A *(p + )] + F A (p^F A (p + ) , 
-Re[G A (p_)G A *(p + )} - F\p_)F\p + ) , 
Re[G A (p_)G A (p + )] - F A (p_)F A (p + ) , 
lm[G A (p^G A *(p + ) = -L A 1 (p;q), 
ReG A (p.)F\ P+ ) + F A (p.)ReG A ( P+ ) = L A (p; q) , 
ImG A (p-)F A ( P+ ) - F A (p_)ImG A ( P+ ) = -L A 2 (p;q) , 



(A6) 



where p + = p + q/2, p_ = p — q/2. 

The matrices representing the Nambu vertex in our approximation for the Hubbard model 
are not all full, that is, several matrix elements vanish. More generally, for coupling functions 
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V PH > A (A;,A;';g) = 



(A7) 



with a factorized momentum dependence and even parity form factors the (direct) particle- 
hole contribution to the vertex has the form 

[2M£ k ,(q) \ 

AUq) P@(q) 2X&(q) 

V 2X&(q) 2X'£(q) 2C A ,(g) J 

and the particle-particle contribution V PP,A has only diagonal elements given by the mag- 
netic coupling function M kk ,(q). However, the matrix elements of the crossed particle-hole 
contribution V PH '' A (A;, k'; q) are all non-zero and generally given by linear combinations of 
several coupling functions. 
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